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We study the asymptotic expansion of the neutral-atom energy as the atomic number Z ^ oo, 
presenting a new method to extract the coefficients from oscillating numerical data. We find that 
recovery of the correct expansion is an exact condition on the Kohn-Sham kinetic energy that is 
important for the accuracy of approximate kinetic energy functionals for atoms, molecules and 
solids, when evaluated on a Kohn-Sham density. For example, this determines the small gradient 
limit of any generalized gradient approximation, and conflicts somewhat with the standard gradient 
expansion. Tests are performed on atoms, molecules, and jellium clusters. We also give a modern, 
highly accurate parametrization of the Thomas-Fermi density of neutral atoms. 

PACS numbers: 



I. INTRODUCTION 

Ground-state Kohn-Sham (KS) density functional the- 
ory (DPT) is a widely-used tool for electronic struc- 
ture calculations of atoms, molecules, and solids 
in which only the density functional for the exchange- 
correlation energy, i^xcMi must be approximated. But 
a direct, orbital- free density functional theory could be 
constructed if only the non- interacting kinetic energy, Ts , 
were known sufficiently accurately as an explicit func- 
tional of the density [4I . Using it would lead automat- 
ically to an electronic structure method that scales lin- 
early with the number of electrons N (with the possible 
exception of the evaluation of the Hartree energy) . Thus 
the KS kinetic energy functional is something of a holy 
grail of density functional purists, and interest in it has 
recently revived 

In this work, we exploit the "unreasonable accuracy" 
of asymptotic expansions [1] , in this case for large neutral 
atoms, to show that there is a very simple exact condi- 
tion that approximations to Tg niust satisfy, if they are 
to attain high accuracy for total energies of matter. By 
matter, we mean all atoms, molecules, and solids that 
consist of electrons in the field of nuclei, attracted by a 
Coulomb potential. The exact condition is the (known) 
asymptotic expansion of Tg/Z^/^ for neutral atoms, in 
powers of Z~^^^. By careful extrapolation from accu- 
rate numerical calculations up to Z ^ 90, we calculate 
the coefficients of this expansion. We find that the usual 
gradient expansion, derived from the slowly-varying gas, 
but applied to essentially exact densities, yields only a 
good approximation to these coefficients. Thus, all new 
approximations should either build in these coefficients, 
or be tested to see how well they approximate them. 
We perform several tests, using atoms, molecules, jel- 
lium surfaces, and jellium spheres, and analyze two ex- 
isting approximations. In Ref. a related method was 



used to derive the gradient coefficient in modern gen- 
eralized gradient approximations (GGA's) for exchange. 
Given this importance of = ^ 00 as a condition on 
functionals, we revisited and improved upon the exist- 
ing parametrizations of the neutral-atom Thomas-Fermi 
(TP) density. The second-half of the paper is devoted to 
testing its accuracy. 



II. THEORY AND ILLUSTRATION 

Por an A^-electron system, the Hamiltonian is 

= T + Kxt + V^c , (1) 

where T is the kinetic energy operator, Vext the ex- 
ternal potential, and V^o the electron-electron interac- 
tion, respectively. The electron density n(r) yields A^ = 
J cPr n{r), where A^ is the particle number. 

To explain asymptotic exactness, we (re)-introduce the 
C-scaled potential [y] (which is further discussed in Ref. 
0]), given by 



(2) 



where Vcxt (r) is the external potential, and the Thomas- 
Fermi expectation value is K^tW = C^^'^Kxt['^]- In this 
(^-scaling scheme, nuclear positions Rq and charges Za 
of molecules are scaled into C^^/'^Ra and (^Za respec- 
tively. In a uniform electric field, £ — > (^^^^S. For neutral 
atoms, scaling C is the same as scaling Z, and this gives 
Schwinger's asymptotic expansion for the total energy of 
neutral atoms [5, [g] , 
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where cq = 0.768745, a = -1/2, C2 ^ 0.269900, and 
Z is the atomic number. This large Z-expansion gives 
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a remarkably good approximation to the Hartree-Fock 
energy of the neutral atoms, with less than a 10% error 
for H, and less than 0.5% error for Ne. By the virial 
theorem for neutral atoms, T = —E, and T ~ Tg to this 
order in the expansion (since the correlation energy is 
roughly ~ Z). Hence, the non-interacting kinetic energy 
has the following asymptotic expansion. 

Ts = Co + ciZ^ + C2 ^5/3 + . . . (4) 

We say that an approximation to the kinetic energy 
functional is asymptotically exact to the p-th degree if 
it can reproduce the exact cq, ci, . . . , Cp. The three dis- 
played terms in Eq. ([3]) constitute the second-order 
asymptotic expansion for the total energy of neutral 
atoms, and we expect that this asymptotic expansion is 
a better starting point for constructing a more accurate 
approximation to the kinetic energy functional than the 
traditional gradient expansion approximation (GEA). 

The leading term in Eq. (|4]) is given exactly by a local 
approximation to Tg (TF theory), but the leading cor- 
rection is due to higher-order quantum effects, and only 
approximately given by the gradient expansion evaluated 
on the exact density. However, these coefficients are vi- 
tal to finding accurate kinetic energies. Since we know 
that cqZ'^/^ becomes exact as N = Z oo, we define 
ATg = Ts—cqZ"^^^ and investigate ATg as a function of Z. 
How accurate is the asymptotic expansion for ATg? In 
Figure m we evaluate Tg for atoms within the optimized 
effective potential (OEP) Q using the exact exchange 
functional and plot the percentage error in ATg, for all 
atoms and the asymptotic series with just two terms. 
The series is incredibly accurate, with only a 13% error 
for N=2 (He), and 14% for iV=l. Thus, any approxima- 
tion that reproduces the correct asymptotic series (up to 
and including the C2 term) is likely to produce a highly 
accurate Tg. 




Z 



FIG. 1: Percentaee error between ci Z'^ + C2 Z^/^ and ATg = 
Tg - coZ^/\ 

To demonstrate the power and the significance of this 
approach, we apply it directly to the first term (where the 



answer is already known, but perhaps not yet fully appre- 
ciated in the DFT community). Using any (all-electron) 
electronic structure code, one calculates the total ener- 
gies of atoms for a sequence running down a column. 
By sticking with a specific column, one reduces the os- 
cillatory contributions across rows, and the alkali-earth 
column yields the most accurate results. By then fitting 
the resulting curve of T^jZ'' 1'^ as a function of Z~^^^ to a 
parabola, one finds cq = 0.7705. Now assume one wishes 
to make a local density approximation (LDA) to Tg, but 
knows nothing about the uniform electron gas. Dimen- 
sional analysis yields [l3| 



r(o)[n] = ^g/, / = 




but does not determine the constant, Ag. A similar fit- 
ting of /, based on the self-consistent densities evaluated 
using the OEP exact exchange functional, gives a leading 
term of 0.2677 Z'^/^, yielding ^g = 2.868. Thus we have 
deduced the local approximation to the non-interacting 
kinetic energy. 

A careful inspection of the above argument reveals 
that the uniform electron gas is never mentioned. As 
'N grows, the wavelength of the majority of the parti- 
cles becomes short relative to the scale on which the po- 
tential is changing, loosely speaking, and semiclassical 
behavior dominates. The local approximation is a uni- 
versal semiclassical result, which is exact for a uniform 
gas simply because that system has a constant poten- 
tial. On the basis of that argument, we know the ex- 
act value is ^g = (3/10)(37r2)2/3 = 2.871, demonstrating 
that (for this case) our result is accurate to about 0.1%. 
This argument tells us that the reliability of the local 
approximation is no indicator of how rapidly the density 
varies. That this argument is correct for neutral atoms 
was carefully proven by Lieb and Simon in 1973 [ll[ and 
later generalized by Lieb to all matter 

The focus of the first part of this paper is on the re- 
maining two known coefficients (ci and C2) and how well 
the GEA performs for them. We evaluate those gradi- 
ent terms by fitting asymptotic series exactly and we 
find that the gradient expansion does well, but is not 
exact. From this information, we develop a modified gra- 
dient expansion approximation that reproduces the exact 
asymptotic coefficients c\ and C2, merely as an illustra- 
tion of the power of asymptotic exactness. We test it on 
a variety of systems, finding the expected behavior. 

In Section |Vl we present a parametrization of the 
TF density which is more accurate than previous 
parametrizations. The TF density has a simple scaling 
with Z and becomes relatively exact and slowly-varying 
for a neutral atom as Z — > cxd, breaking down only near 
the nucleus and in the tail. We compare various quanti- 
ties of our parametrization with exact values and earlier 
parametrizations, and analyze the properties of the TF 
density. 
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III. LARGE Z METHODOLOGY 

We begin with a careful methodology for extracting 
the asymptotic behavior from highly accurate numerical 
calculations. Fully numerical DFT calculations were per- 
formed using the OPMKS code [Tl to calculate the total 
energies of neutral atoms using the OEP exact exchange 
functional. The spin-density functional version of Ts has 
been used for all systems [l3l l. 




FIG. 2: Difference between T^jZ^I'-^ and co -I- ci Z'^'^ + 
C2 as a function of with exact asymptotic co- 

efficients. 



To attain maximum accuracy for ci and C2, we need to 
suppress the oscillations which come from the next term, 
C'iZ'^l^ . Consider first the OEP results. We investigate 
the differences between T^^^ IZ'^I'^ and Cq + ci Z^^/^ 
C2 with exact asymptotic coefficients in Figure [21 

We extract 6 data points (Z=24 (Cr), 25 (Mn), 30 (Zn), 
31 (Ga), 61 (Pm), and 74 (W)) which have the smallest 
differences, i.e., nearest to where the curve crosses the 
horizontal axis. We then make a least-squares fit with a 
parabolic form in Z~^^^ . ignoring the oscillation term. 



^ = 0.768745 ^ cxZ-^l"^ + c^Z 



-2/3 



(6) 



Effectively, we solve two linear equations for ci and C2. 
We explicitly include the exact cq — 0.768745, since we 
don't have enough data points to extract cq accurately, 
especially in the region Z~^l'^ < 0.2. It is important to 
control the behavior of the fitting line at Z ^ 00. This 
fitting yields a good estimate of ci = —0.5000 and C2 — 
0.2702, with error less than 1%. This demonstrates the 
accuracy of our method for ci and C2 (by construction). 

We repeat the same procedure to extract ci and C2 
coefficients of TF and second- and fourth-order GEA's 
which are given by 



and @,[Tl[il: 

2nGEA4 _ rpTF rp{2) y(4) 



(8) 



These gradient corrections to the local approximation are 
given by 



and 

y(4) 



81 



g2(r)-|g(r).2(r)-ff!M 



(9) 



where T"^^(r), s(r), and q{r) are defined as 



s{r) 



IVnfr) 



2fcp(r)n(r) 



^ V^n(r) 
' 4fc|(r)n(r) 



(11) 



(12) 



(13) 



and fcp(r) {3TT^n{r))^/^ . 

We have also applied this procedure to both T^^^ and 
T^^\ Since the asymptotic expansions of these energies 
begin at Z^, we extract only a ci and a C2 for each using 
the following equations: 



2iGEA2 _ rpTF 

ZV^ 

zv'^ 



= -f AC2Z-2/3 . (14) 



T 



GEA2 



T 



TF 



(7) 



These results are also included in Table HI and are of 
course consistent with our results from Eq. 



IV. RESULTS AND INTERPRETATION 

To understand the meaning of the above results, begin 
with the values of ci. We have combined the results of 
the T(2) and T^'*) fits with that of the T^f fit to produce 
the asymptotic coefficients of T^^^^ and T^^^^. We 
check that these combinations produce the same coeffi- 
cients in TableUwhich are found from the direct fitting of 

2nGEA2 yGEA4 ^gj^^g ^ rpj^^ ^^^^^ ^^-^^^ 

is —1/2. We see that the local approximation (TF) gives 
a good estimate, —0.66. Then the second-order gradient 
expansion yields —0.54, reducing the error by a factor 
of 5. Finally, the fourth-order gradient expansion yields 
—0.52, a further improvement, yielding only a 4% error 
in its approximation to the Scott correction (l6j . 

For C2, the gradient expansion is less useful. The ex- 
act result is 0.27, while the TF approximation overes- 
timates this as 0.39. The GEA2 result is only slightly 
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Cl 


C2 


Exact 


-0.5000 


0.2699 


yOEP 


-0.5000 


0.2702 


rpTF 


-0.6608 


0.3854 


y(2) 


0.1246 


-0.0494 


y(4) 


0.0162 


0.0071 


yGEA2 


-0.5362 


0.3360 


J.GEA4 


-0.5200 


0.3431 


yGGAa 


-0.5080 


0.2918 


j^LmGGAa 


-0.5089 


0.3174 



"See section ITVl 

TABLE I: The coefficients in the asymptotic expansion of the 
exact kinetic energy and various local and semilocal function- 
als. The fit was made to Z=24 (Cr), 25 (Mn), 30 (Zn), 31 
(Ga), 61 (Pm), and 74 (W). The functionals of the last two 
rows are defined in section HVl 



can be improperly negative for rapidly-varying densities, 
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reduced (0.34), and the fourth-order correction has the 
wrong sign. 

To understand how important these results can be, we 
consider how exchange and correlation functionals are 
constructed. Often, such constructions begin from the 
GEA, which is then generalized to include (in an approx- 
imate way) all powers of a given gradient. For slowly 
varying densities, it is considered desirable to recover the 
GEA result. But we have seen here how this conflicts 
with the asymptotic expansion, and in Ref [5], it was 
shown how the asymptotic expansion is more significant 
to energies of real materials, and how successful GGA's 
for atoms and molecules well- approximate the large-Z 
asymptotic result, not the slowly- varying gas. 

Atoms: To illustrate this point, we construct here a 
trivial modified gradient expansion, MGEA2, designed 
to have the correct asymptotic coefficients, in so far as is 
possible. Thus 

j,MGEA2^yTF^^ 290r(2) (15) 

The enhancement coefficient has been chosen to make 
^MGEA2 = _i/2 exactly. In Table Ull we list the results 
of several different approximations for the alkali-earth 
atoms. Because the GEA2 error passes through around 
Z=8, its errors are artificially low. 

We can repeat this exercise for the fourth order, re- 
quiring both Cl and C2 be exact. Now we find: 

T^^^^^[n] = T^^[n] + 1.789 T^^^[n]- 3. MlT^^^n] (16) 

i.e., strongly modified gradient coefficients. This is some- 
what arbitrary, as there are several terms in T^^-* , and 
there's no real reason to keep their ratios the same as in 
GEA (Eq. (Uni)). However, the results of Table |TT] and 
Figure [3] speak for themselves. The resulting functional 
is better than either GEA for all the alkali-earths. Of 
course, the exact is positive for any density, as are the 
terms T'^^ , T^^) and T^^^ of the GEA. Eq. dllD however 



FIG. 3: Percentage errors for atoms (from Z = 1 to Z — 92) 
using various approximations. 



Molecules: The improvement in total kinetic energies 
is not just confined to atoms. Also, for non-interacting 
kinetic energies of molecules, using the data in Ref. p7j . 
Eq. gives better average of the absolute errors in 

hartree (0.6) than T^f (9.4), T'^^^^ (0.9), and TC^eA4 
(0.8), shown in Table ITTTl Of greater importance are en- 
ergy differences. For atomization kinetic energies, also 
using the data in Ref. [17], T"^^ gives the best averaged 
absolute error (0.25), which is worsened by gradient cor- 
rections. Since the GEA does not have the right quantum 
corrections from the edges, turning points and Coulomb 
cores , GEA does not improve on the atomization pro- 
cess. However, the TF kinetic energy functional is always 
the dominant term. So, TF gives very good results on 
the atomization kinetic energies. But the error (0.29) 
of Eq. (Uni) is smaller than that of T'^^^^ (0.36) and 
yGEA4 (0.44). In either case, Eq. ^ works better for 
atoms and molecules than the fourth-order gradient ex- 
pansion. Thus, requiring asymptotic exactness is a useful 
and powerful constraint in functional design. 

Jellium surfaces: We test this MGEA4 functional 
for jellium surface kinetic energies. As shown in Table 
IIVI the T*^"*^ term in T^^^^ improves the jellium surface 
kinetic energy in comparison to the results of rGEA2 
but Eq. (jl6p worsens the jellium surface kinetic energies 
due to the strongly modified coefficient of T'^^\ This is 
a confirmation of our general approach. By building in 
the correct aysmptotic behavior for atoms, including the 
Scott correction coming from the Is region, we worsen 
energetics for systems without this feature. 

Jellium spheres: We also investigate the kinetic en- 
ergies of neutral jellium spheres (with KS densities using 
LDA exchange-correlation and with rg = 3.9) from Ref. 
. The analysis of the results is based upon the liquid 
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Atom 


Z 


yOEP 


rpTF 


%err 


j,GEA2 


%err 


yMGEA2 


%err 


rpGEAl 


%err 


yMGEA4 


%err 


Be 


4 


14.5724 


13.1290 


-10 


14.6471 


0.5 


15.0880 


3.5 


14.9854 


2.8 


14.5453 


-0.2 


Mg 


12 


199.612 


184.002 


-8 


198.735 


-0.4 


203.014 


1.7 


201.452 


0.9 


199.924 


0.2 


Ca 


20 


676.752 


630.064 


-7 


672.740 


-0.6 


685.136 


1.2 


680.286 


0.5 


677.433 


0.1 


Sr 


38 


3131.53 


2951.89 


-6 


3110.44 


-0.7 


3156.50 


0.8 


3136.76 


0.2 


3134.48 


0.09 


Ba 


56 


7883.53 


7478.27 


-5 


7829.36 


-0.7 


7931.34 


0.6 


7886.19 


0.03 


7888.14 


0.06 


Ra 


88 


23094.3 


22065.8 


-4 


22945.9 


-0.6 


23201.5 


0.5 


23083.9 


-0.05 


23110.5 


0.07 



TABLE II; KS kinetic energy (T) in hartrees and various approximations for alkali-earth atoms. 



Atom 




rpTFa 


j,GEA2a 


rpGEAia 


rpMGEA4 


n 


u.ouu 


n 0/1/1 


U.Ull 


u.uoz 


-U.UZD 


JD 


OA KAQ 




-u.uuo 


n A7(\ 


n 1 77 

-U. 1 t t 




■^7 71 A 


-0. / 01 


n 1 


n finn 


n 998 


NT 




/I 00*^ 


n nQ7 




n 078 




74 8fi7 


-u. y yu 






n AQ7 


r 




-y.uyo 


-u.yoo 


u.Doy 


-u.Duy 


112 


1.151 


-0.142 


-0.014 


0.033 


-0.094 


HP 


100.169 


-9.016 


-0.920 


0.639 


-0.520 


H2O 


76.171 


-7.074 


-0.692 


0.565 


-0.484 


CH4 


40.317 


-3.773 


-0.140 


0.619 


-0.189 


NH3 


56.326 


-5.292 


-0.400 


0.587 


-0.331 


BF3 


323.678 


-29.052 


-2.641 


2.454 


-1.370 


CN 


92.573 


-8.940 


-0.687 


0.978 


-0.570 


CO 


112.877 


-10.694 


-0.911 


1.036 


-0.670 


F2 


199.023 


-18.367 


-2.201 


0.925 


-1.451 


HCN 


92.982 


-8.925 


-0.658 


1.008 


-0.534 


N2 


109.013 


-10.487 


-0.916 


0.999 


-0.719 


NO 


129.563 


-12.342 


-1.240 


0.962 


0.279 


O2 


149.834 


-14.186 


-1.527 


0.965 


-1.110 


O3 


224.697 


-21.636 


-2.699 


1.028 


-2.071 


MAE' 




9.364 


0.872 


0.812 


0.600 



"Ref. [T3l 

''Mean absolute error 

TABLE III: Exact non-interacting kinetic energy (in hartrees) 
for molecules, and errors in approximations. All values are 
evaluated on the converged KS orbitals and densities obtained 
with B88-PW91 functionals, and the MGEA4 kinetic energies 
are evaluated using the TP and the GEA data from Ref. [17| . 

drop model of Refs. [li,[23|. We write 

Ts{rs,N) = ^7ri?V-"f (rs) + ^nR^a, + 27rRjf {r^, N), 

(17) 

where R is the radius of the sphere of uniform positive 
background. Since we know the bulk (uniform) kinetic 
energy density, r^"'^, and the surface kinetic energy CTs 
for a given functional, we can extract 7|*^(rs,A'^) from 
this equation, and 

hm jf{r„N)=j,{r,) (18) 

N^oo 

is the curvature energy of jellium. We calculate 



rs 


Exact 


rpTF 


21GEA2 


21GEA4 


yMGEA2a 


yMGEA4 6 


j.LmGGA 


2 


-5492.7 


11 


2.5 


1.1 


-0.9 


0.73 


1.3 


4 


-139.9 


54 


22 


11 


12 


36 


15 


6 


-3.4 


660 


330 


180 


238 


675 


280 



"See Eq. ^ 
'See Eq. ^ 



TABLE IV: Exact jellium surface kinetic energies [erg/cm^) 
and % error, which is ((j|^^ — (t™)/^™, of each approximation. 

^f{rs,N) using the TF, GEA, MGEA, and a Laplacian- 
level meta-GGA (LmGGA) of Ref. [18(1, which is ex- 
plained further in the following subsection. From Ta- 
ble |Vl we observe that: (i) Gradient corrections in GEA 



N 


Exact 


yGEA2 


yGEA4 


yMGEA2a 


rpMGEAlb 


jiLmGGA 


2 


-1.8 


1.1 


2.4 


1.5 


-2.8 


1.9 


8 


-1.9 


1.0 


2.1 


1.3 


-2.3 


-5.1 


18 


-0.5 


1.2 


2.0 


1.6 


-0.7 


-6.4 


58 


-0.8 


1.3 


2.2 


1.7 


-1.1 


-3.2 


92 


-1.7 


1.2 


2.0 


1.5 


-1.0 


-1.9 


254 


-0.5 


1.4 


2.3 


1.8 


-0.9 





"See Eq. lO 
'See Eq. l(T6)l 



TABLE V: 10''x(7f (rs,iV)-7j^(rs,Af)) in atomic units vs. 
N = Z for neutral jellium spheres with rs = 3.93 with various 
functionals. As A'^ = Z — > 00, 71"' tends to the curvature 
kinetic energy of jellium, 73. 

worsen 7^^. (ii) The LmGGA of Ref. [3 is even worse 
than rGEA4^ ^ (which has the right Cq and 

ci) is not so good, but better than T^^^''. (iv) Eq. HH) 
(which has the right cq, ci, and C2) gives good results. 

Existing approximations: We suggest that the 
large-Z asymptotic expansion is a necessary condition 
that an accurate kinetic energy functional should satisfy, 
but is not sufficient. We show this by testing two kinds of 
semilocal approximations (GGA and meta-GGA) to the 
kinetic energy functionals. 

Recently, Tran and Wesolowski [2l| constructed a 
GGA-type kinetic energy functional using the conjoint- 
ness conjecture. They found the enhancement factor 
by minimizing mean absolute errors of kinetic energies 
for closed-shell atoms. We evaluate the kinetic energies 
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of atoms using this functional (tGGA-j extract the 
asymptotic coefficients shown in Table HI This gives a 
good ci coefficient, with C2 close to the exact value, and 
so is much more accurate than the GEA's. 

Perdew and Constantin [18] constructed a LmGGA for 
the positive kinetic energy density r that satisfies the lo- 
cal bound T > Tw, where tw is the von Weizsacker ki- 
netic energy density, and tends to tw as r ^ in an 
atom. It recovers the fourth-order gradient expansion in 
the slowly- varying limit. We calculate the asymptotic 
coefficients shown in Table U for this functional. These 
values are better than those of T'^^^"'. The good ci from 
2iGEA4 appears somewhat fortuitous, since there is noth- 
ing about a slowly-varying density that is relevant to a 
cusp in the density. The good Scott correction ci from 
the LmGGA comes from correct physics: LmGGA recov- 
ers the von Weizsacker kinetic energy density in the Is 
cusp, without the spurious but integrable divergences of 
the integrand of T'^^^^. 

We finish by discussing other columns of the periodic 
table. We have also performed all these calculations on 
the noble gases. In fact, from studies of the asymptotic 
series [2^ . it is known that the shell-structure occurs in 
the next order, Z'^^^, and that the noble gases are furthest 
from the asymptotic curves. But Table IVII shows our 
functionals work almost as well for the noble gas series. 



V. MODERN PARAMETRIZATION OF 
THOMAS-FERMI DENSITY 

Our asymptotic expansion study gives new reasons for 
studying large Z atoms. Our approximate functionals 
were tested on highly accurate densities, but ultimately, 
self-consistency is an important and more-demanding 
test. Any approximate functional yields an approxi- 
mate density via the Euler equation. In this section, 
we present a new, modern parametrization of the neutral 
atom TF density, which is more accurate than earlier 
versions [l^, [ISl ■ 

The TF density of a neutral atom can be written as 



n{r) 



3/2 



(19) 



where a = (l/2)(37r/4)2/3 and x = Z^/^r/a, and the 
dimensionless TF differential equation is 



$3(a;) 



$(a;) > 0, 



(20) 



dx'^ V X 
which satisfies the following initial conditions: 

$(0) = 1, $'(0) = -S, 5 = 1.5880710226. (21) 

We construct a model for <& which recovers the first eight 
terms of the small- a; expansion and the leading term of 
the asymptotic expansion at large- a; {^{x) lAA/x^, 
as X — > oo). Following Tal and Levy [25|, we use y = 



y/x as the variable, because of the singularity of the TF 
equation. Our parametrization is 




y=x 



FIG. 4: The exact numerical ${y) and parametrized <I>(y) can 
not be distinguished. 




10 12 14 



y = x 



FIG. 5: Errors in the model, relative to numerical integra- 
tions. 




$™°'*(y) 



where and /3i are coefficients given in the Table I VIII 
The values of ai arc fixed by the small y-expansion, while 
those of f3i are found by minimization of the weighted 
sum of squared residuals, x^j for < y < 10. The was 
minimized using the Levenberg-Marquardt method f23|. 
This method is for fitting when the model depends non- 
linearly on the set of unknown parameters. 1000 points 
were used, equally spaced between y = and y = 10. We 
plot the numerically exact ^{y) and our model in Figure 
m and the differences between them in Figure [51 These 
graphs illustrate the accuracy of our parametrization. 
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Atom 


Z 


yOEP 


rpTF 


%err 


yGEA2 


%err 


yMGEA2 


%err 


j,GEA4 


%err 


yMGEA4 


%err 


He 


2 


2.86168 


2.56051 


-11 


2.87847 


0.6 


2.97083 


3.8 


2.96236 


3.5 


2.80717 


-1.9 


Ne 


10 


128.545 


117.761 


-8 


127.829 


-0.6 


130.753 


1.7 


129.737 


0.9 


128.447 


-0.08 


Ar 


18 


526.812 


489.955 


-7 


524.224 


-0.5 


534.178 


1.4 


530.341 


0.7 


527.772 


0.2 


Kr 


36 


2752.04 


2591.20 


-6 


2733.07 


-0.7 


2774.27 


0.8 


2756.72 


0.2 


2754.17 


0.08 


Xe 


54 


7232.12 


6857.94 


-5 


7183.78 


-0.7 


7278.42 


0.6 


7236.65 


0.06 


7237.85 


0.08 


Rn 


86 


21866.7 


20885.7 


-4 


21725.4 


-0.6 


21969.3 


0.5 


21857.2 


-0.04 


21881.7 


0.07 



TABLE VI: KS kinetic energy (T) in hartrees and various approximations for noble atoms. 



a2 


-B 


/3i 


-0.0144050081 


as 


4/3 




0.0231427314 


as 


-2B/5 


Pa 


-0.00617782965 


ae 


1/3 


Pa 


0.0103191718 


ar 


3BV70 


P5 


-0.000154797772 


as 


-2B/15 






ag 


2/27 + 







TABLE VII: The values of Pi are found by fitting Ea. (p2)l 
to the exact solution, and those of ai are the parameters of 
small-y expansion [23 • B is given by 1.5880710226. 



To compare the quality of the various parametriza- 
tions, we calculate the p-th moment of the j-th power 
of ^{x)/x: 



(p) 



dxx^ 



(28) 



Many quantities of interest can be expressed in terms of 
these moments: 

1) Particle number: To ensure we 
require 



In Table IVIIII we calculate several moments using 
our model and existing models that were proposed by 
Gross and Dreizler [23| and Latter The Latter 

parametrization is 

$^(a;) = 1/(1 0.02747x^/2 -h 1.243a; - 0.1486x^/2 
-h0.2303x2 + 0.007298x^/2 
+0.006944x3) ^ (23) 

and the Gross-Dreizlcr model (which correctly removes 
the ^/x term) is: 



q>GD{x) = 1/(H- 1.4712X- 0.4973x^/2 
0.3875x2 -f 0.002102x3) . 



(24) 



Lastly, we introduce an extremely simple model that we 
have found useful for pedagogical purposes (even when 
N differs from Z). We write 



ped 



(r) 



N 



1 



27r3/2i?3/2 j.3/2 



-r/R 



R 



Z - (5N 



(25) 



where a = (9/5^/5)(\/37^/4)l/3 and /3 = 1/2 - I/tt 
have been found from integration of the TF kinetic and 
Hartree energies, respectively, and R minimizes the TF 
total energy. For N = Z, this yields: 



(26) 



ci>P-d., -2ail-P)./3a ( I ^ i] 

^ ^ ^ ' ' 6V3 V2 ^/ 

This crude approximation does not satisfy the correct 
initial conditions of Eq. pTjl : 



^P'"'{0) ^ 7 = 0.880361(^1), 



125(2- 



648(47r5)i/3 



-0.48 -1.59). (27) 



Mi% = 1 



(29) 



2) TF kinetic energy: The TF kinetic energy is cq^^/^, 
which implies 



M, 



(2) 



-B . 



5/2 J 

3) The Hartree energy is U ^ ^ J J d^r d^r' ^4^^ 



(30) 



^M^/^Z^/3, which imphes 



(31) 



energy is defined as T4xt = 
iM^J^Z''/3 for the exact TF den- 
3D. 



4) The external 
— J d'^r Z n{r)/r = - 

sity, which also implies Eq, 

5) The local density approximation (LDA) ex- 
change energy is defined as E^^^ — J dr n^l'^ (r) 
where Ax = -(3/4)(3/7r)i/3, so for TF, = 
ylx(47ra3)("i/3)7v/f ^^5/3, which implies 



uf^ = 0.615434679 . 



(32) 



(2) 

This M2 is evaluated on the exact TF density which 
we calculate numerically. LDA exchange suffices 1^, |5| 
for asymptotic exactness to the order displayed in Eqs. 
^ and ([U; for a numerical study, see Ref. [13] ■ Table 
EnT] shows that our modern parametrization is far more 
accurate than existing models by all measures, and that 
our simple pedagogical model is roughly correct for many 
features. 
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TABLE VIII: Various moments calculated with our model and with the models of Ref. [2^, [2J. Here M'^^^ is given by 



moment 


our model 


% error 


Gross and Dreizler [23] 


% error 


Latter [24] 


% error 




% error 


exact 


-'"3/2 


0.999857885 


-0.01 


1.008 


0.8 


0.999 


-0.04 


1 





1 


^"^5/2 


1.13426462 


-0.006 


1.1299 


-0.4 


1.137 


0.2 


1.11 


-2 


5B/7 




0.615438208 


0.001 


0.6129 


-0.4 


0.616 


0.02 


0.72 


16 


0.615434679" 




1.58799857 


-0.005 


1.5844 


-0.2 


1.589 


0.07 


1.62 


2 


B 



"Numerical result from the TF differential equation. 



Finally, we make some comparisons with densities of 
real atoms to illustrate those features of real atoms that 
are captured by TF. The radial density, s(r) (Eq. (fT^). 
and q{r) (Eq. p3)) ) are given by 



where f{x) — \/x'^'^/'^{x), 



(33) 



KO-^^, «i = (9/2^)^/^2, (34) 



and 



q{r) 



al {52(0;) +2x2$(x)$"(x)} 



3Z2/3 



(35) 



where g{x) is defined as ^(x) — a;$'(x). The gradient 
relative to the screening length is 

t{v) = 2^^^!^ , where k, (r) = V4fcp(r)/^ , (36) 



and here 



a2\g{x)\ 
(a;3$5(x))i/4 



02 



35^6^ 

28/3. A 



= 0.6124. 



(37) 



We also show large- and small- a; limit behaviors of var- 
ious quantities using <^{x) \AA/x^ as x — > 00 and 
<^{x) I - Bx^ as a; ^ 0. 



Z2 1 


a;-^0 






Z4/3 ^ 




y/x 




a 




ai 1 


x^O 






al 1 


x^O 


3Z2/3 x 




a2 


x^O 


3.3/4 





n{r) 

s{r) 

q{r) 
tir) 



432^2 



144^4/3 

aix 
3Zi73 ' 
balx"^ 

2a2 
V3- 



(38) 

(39) 
(40) 

(41) 
(42) 



We plot the Z-scaled exact (self-consistent densities 
with OEP exact exchange functional) and TF radial den- 
sities of Ba {Z = 56) and Ra {Z = 88) in Figure El Al- 
though the shell structure is missing, and the decay at 




FIG. 6: Plot of the scaled radial densities of Ba and Ra using 
Eq.(l33} and SCF densities with OEP exact exchange. TF 
scaled densities of Ba and Ra are on top of each other. 




FIG. 7: Plot of the scaled reduced density gradient a(r) (rel- 
ative to the local Fermi wavelength) vs. Z^^^r. 



a large distance is wrong, the overall shape of the TF 
density is relatively correct. 

In Figures [3 [3 and [H we plot the scaled s(r), q[r), 
and t{r) using the exact and the TF densities of Ba and 
Ra. In particular, t(r) measures how fast the density 
changes on the scale of the TF screening length, and its 



9 



magnitude does not vary with Z in TF theory. From 
these figures, we see that sir), qir) and i(r) of the TF 
density diverge near the nucleus, since the TF density 
does not satisfy Kato's cusp condition. 




FIG. 8; Plot of the scaled reduced Laplacian q(r) (relative to 
the local Fermi wavelength) vs. Z^^^r. 



prove useful in the never-ending search for improved den- 
sity functionals. 




FIG. 9: Plot of the reduced density gradient t{r) (relative to 
the local screening length) vs. Z^^^r. As r ^ oo, the TF 
t 0.7071. 



When N = Z ^ oo for a reahstic density, s(r) is small 
except in the density tail (s ^ Z~^l^ over most of the 
density), and g(r) is small except in the tail and Is core 
regions (g ~ Z~'^^^ over most of the density). This is 
why gradient expansions for the kinetic and exchange 
energies, applied to realistic densities, work as well as 
they do in this limit. The kinetic and exchange energies 
have only one characteristic length scale, the local Fermi 
wavelength, but the correlation energy also has a different 
one, the local screening length. Since t(r) is not and does 
not become small in this limit, gradient expansions do not 
work well at all for the correlation energies of atoms Q. 
The standard of "smallness" for s and g, and the more 
severe standard of smallness for are explained in Refs. 
and m. 

Finally we evaluate T^") + T^^) on the TF density. We 
find the correct cq in the Z ^ oo expansion from T'^''^ 
but ci vanishes, due to the absence of a proper nuclear 
cusp, and C2 diverges because T^^^ diverges at its lower 
limit of integration. 

VI. SUMMARY 

We have shown the importance of the large-A^ limit 
for density functional construction of the kinetic energy 
(with the functional evaluated on a Kohn-Sham density), 
and also provided a modern, highly accurate parameteri- 
zation of the neutral-atom TF density. Our results should 



For atoms and molecules, the large- A'^ limit seems more 
important than the slowly- varying limit. On the ladder 
[291] of density-functional approximations, there are three 
rungs of semilocal approximations (followed by higher 
rungs of fully nonlocal ones). The LDA uses only the 
local density, the GGA uses also the density gradient, 
and the meta-GGA uses in addition the orbital kinetic 
energy density or the Laplacian of the density. For the 
exchange-correlation energy, the GGA rung cannot [^.[281 
simultaneously describe the slowly-varying limit and the 
iV = Z — > 00 limit for an atom, and we have found here 
that the same is true (but less severely by percent error 
of a given energy component) for the kinetic energy. This 
follows because, as iV = Z ^ 00, the reduced gradient 
sir) of Eq. (|12p becomes small over the energetically im- 
portant regions of the atom, as can be inferred from Fig. 
[3 so that a GGA reduces to its own second-order gra- 
dient expansion even in regions where a meta-GGA does 
not [51] (e.g., near a nucleus, where q[r) diverges but s(r) 
does not, as shown in Figs. [7] and [8]). For the kinetic as 
for the exchange-correlation energy, meta-GGA's jlsl can 
recover both the slowly-varying and large-Z limits; it re- 
mains to be seen how well fully nonlocal approximations 
[30I [sH can do this. 

We thank Eberhard Engcl for the use of his atomic 
OPMKS code, and NSF (CHE-0355405 and DMR- 
0501588) and the Korea Science and Engineering Foun- 
dation Grant (No. C00063), for support. 
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